// Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.

#include <vector>
#include "GaussIntegration.h"
#include "GaussLegendre1D.h"

IntPt GQT1[1] = {
  {{0.333333333333333, 0.333333333333333, 0.}, 0.500000000000000}};
IntPt GQT2[3] = {
  {{0.166666666666667, 0.166666666666667, 0.}, 0.166666666666667},
  {{0.666666666666667, 0.166666666666667, 0.}, 0.166666666666667},
  {{0.166666666666667, 0.666666666666667, 0.}, 0.166666666666667}};

IntPt GQT3[4] = {
  {{0.333333333333333, 0.333333333333333, 0.}, -0.281250000000000},
  {{0.600000000000000, 0.200000000000000, 0.}, +0.260416666666667},
  {{0.200000000000000, 0.600000000000000, 0.}, +0.260416666666667},
  {{0.200000000000000, 0.200000000000000, 0.}, +0.260416666666667}};

IntPt GQT4[6] = {
  {{0.816847572980459, 0.091576213509771, 0.}, 0.054975871827661},
  {{0.091576213509771, 0.816847572980459, 0.}, 0.054975871827661},
  {{0.091576213509771, 0.091576213509771, 0.}, 0.054975871827661},
  {{0.108103018168070, 0.445948490915965, 0.}, 0.111690794839005},
  {{0.445948490915965, 0.108103018168070, 0.}, 0.111690794839005},
  {{0.445948490915965, 0.445948490915965, 0.}, 0.111690794839005}};

IntPt GQT5[7] = {
  {{0.333333333333333, 0.333333333333333, 0.}, 0.112500000000000},
  {{0.797426985353087, 0.101286507323456, 0.}, 0.062969590272414},
  {{0.101286507323456, 0.797426985353087, 0.}, 0.062969590272414},
  {{0.101286507323456, 0.101286507323456, 0.}, 0.062969590272414},
  {{0.470142064105115, 0.059715871789770, 0.}, 0.066197076394253},
  {{0.059715871789770, 0.470142064105115, 0.}, 0.066197076394253},
  {{0.470142064105115, 0.470142064105115, 0.}, 0.066197076394253}};

IntPt GQT6[12] = {
  {{0.873821971016996, 0.063089014491502, 0.}, 0.025422453185104},
  {{0.063089014491502, 0.873821971016996, 0.}, 0.025422453185104},
  {{0.063089014491502, 0.063089014491502, 0.}, 0.025422453185104},
  {{0.501426509658179, 0.249286745170910, 0.}, 0.058393137863189},
  {{0.249286745170910, 0.501426509658179, 0.}, 0.058393137863189},
  {{0.249286745170910, 0.249286745170910, 0.}, 0.058393137863189},
  {{0.636502499121399, 0.310352451033785, 0.}, 0.041425537809187},
  {{0.310352451033785, 0.636502499121399, 0.}, 0.041425537809187},
  {{0.636502499121399, 0.053145049844816, 0.}, 0.041425537809187},
  {{0.310352451033785, 0.053145049844816, 0.}, 0.041425537809187},
  {{0.053145049844816, 0.310352451033785, 0.}, 0.041425537809187},
  {{0.053145049844816, 0.636502499121399, 0.}, 0.041425537809187}};

IntPt GQT7[13] = {
  {{0.333333333333333, 0.333333333333333, 0.}, -0.074785022233841},
  {{0.479308067841920, 0.260345966079040, 0.}, +0.087807628716604},
  {{0.260345966079040, 0.479308067841920, 0.}, +0.087807628716604},
  {{0.260345966079040, 0.260345966079040, 0.}, +0.087807628716604},
  {{0.869739794195568, 0.065130102902216, 0.}, +0.026673617804419},
  {{0.065130102902216, 0.869739794195568, 0.}, +0.026673617804419},
  {{0.065130102902216, 0.065130102902216, 0.}, +0.026673617804419},
  {{0.048690315425316, 0.312865496004874, 0.}, +0.038556880445128},
  {{0.312865496004874, 0.048690315425316, 0.}, +0.038556880445128},
  {{0.638444188569810, 0.048690315425316, 0.}, +0.038556880445128},
  {{0.048690315425316, 0.638444188569810, 0.}, +0.038556880445128},
  {{0.312865496004874, 0.638444188569810, 0.}, +0.038556880445128},
  {{0.638444188569810, 0.312865496004874, 0.}, +0.038556880445128}

};

IntPt GQT8[16] = {
  {{0.333333333333333, 0.333333333333333, 0.}, 0.072157803838894},
  {{0.081414823414554, 0.459292588292723, 0.}, 0.047545817133643},
  {{0.459292588292723, 0.081414823414554, 0.}, 0.047545817133643},
  {{0.459292588292723, 0.459292588292723, 0.}, 0.047545817133643},
  {{0.658861384496480, 0.170569307751760, 0.}, 0.051608685267359},
  {{0.170569307751760, 0.658861384496480, 0.}, 0.051608685267359},
  {{0.170569307751760, 0.170569307751760, 0.}, 0.051608685267359},
  {{0.898905543365938, 0.050547228317031, 0.}, 0.016229248811599},
  {{0.050547228317031, 0.898905543365938, 0.}, 0.016229248811599},
  {{0.050547228317031, 0.050547228317031, 0.}, 0.016229248811599},
  {{0.008394777409958, 0.728492392955404, 0.}, 0.013615157087217},
  {{0.728492392955404, 0.008394777409958, 0.}, 0.013615157087217},
  {{0.263112829634638, 0.008394777409958, 0.}, 0.013615157087217},
  {{0.008394777409958, 0.263112829634638, 0.}, 0.013615157087217},
  {{0.263112829634638, 0.728492392955404, 0.}, 0.013615157087217},
  {{0.728492392955404, 0.263112829634638, 0.}, 0.013615157087217}};

IntPt *GQT[9] = {GQT1, GQT1, GQT2, GQT3, GQT4, GQT5, GQT6, GQT7, GQT8};
int GQTnPt[9] = {1, 1, 3, 4, 6, 7, 12, 13, 16};

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 1 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP1Solin[1] = {
  {{0.333333333333333, 0.333333333333333, 0.}, 0.500000000000000}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 2 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP2Solin[3] = {
  {{0.166666666666667, 0.166666666666667, 0.}, 0.166666666666667},
  {{0.166666666666667, 0.666666666666667, 0.}, 0.166666666666667},
  {{0.666666666666667, 0.166666666666667, 0.}, 0.166666666666667}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 3 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP3Solin[4] = {
  {{0.333333333333333, 0.333333333333333, 0.}, -0.281250000000000},
  {{0.200000000000000, 0.200000000000000, 0.}, 0.260416666666667},
  {{0.200000000000000, 0.600000000000000, 0.}, 0.260416666666667},
  {{0.600000000000000, 0.200000000000000, 0.}, 0.260416666666667}};
// 1 negative weight, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 4 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP4Solin[6] = {
  {{0.445948490915965, 0.445948490915965, 0.}, 0.111690794839005},
  {{0.445948490915965, 0.108103018168070, 0.}, 0.111690794839005},
  {{0.108103018168070, 0.445948490915965, 0.}, 0.111690794839005},
  {{0.091576213509771, 0.091576213509771, 0.}, 0.054975871827661},
  {{0.091576213509771, 0.816847572980459, 0.}, 0.054975871827661},
  {{0.816847572980459, 0.091576213509771, 0.}, 0.054975871827661}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 5 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP5Solin[7] = {
  {{0.333333333333333, 0.333333333333333, 0.}, 0.112500000000000},
  {{0.470142064105115, 0.470142064105115, 0.}, 0.066197076394253},
  {{0.470142064105115, 0.059715871789770, 0.}, 0.066197076394253},
  {{0.059715871789770, 0.470142064105115, 0.}, 0.066197076394253},
  {{0.101286507323456, 0.101286507323456, 0.}, 0.062969590272414},
  {{0.101286507323456, 0.797426985353087, 0.}, 0.062969590272414},
  {{0.797426985353087, 0.101286507323456, 0.}, 0.062969590272414}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 6 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP6Solin[12] = {
  {{0.249286745170910, 0.249286745170910, 0.}, 0.058393137863189},
  {{0.249286745170910, 0.501426509658179, 0.}, 0.058393137863189},
  {{0.501426509658179, 0.249286745170910, 0.}, 0.058393137863189},
  {{0.063089014491502, 0.063089014491502, 0.}, 0.025422453185104},
  {{0.063089014491502, 0.873821971016996, 0.}, 0.025422453185104},
  {{0.873821971016996, 0.063089014491502, 0.}, 0.025422453185104},
  {{0.310352451033785, 0.636502499121399, 0.}, 0.041425537809187},
  {{0.636502499121399, 0.053145049844816, 0.}, 0.041425537809187},
  {{0.053145049844816, 0.310352451033785, 0.}, 0.041425537809187},
  {{0.310352451033785, 0.053145049844816, 0.}, 0.041425537809187},
  {{0.636502499121399, 0.310352451033785, 0.}, 0.041425537809187},
  {{0.053145049844816, 0.636502499121399, 0.}, 0.041425537809187}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 7 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP7Solin[13] = {
  {{0.333333333333333, 0.333333333333333, 0.}, -0.074785022233841},
  {{0.260345966079040, 0.260345966079040, 0.}, 0.087807628716604},
  {{0.260345966079040, 0.479308067841920, 0.}, 0.087807628716604},
  {{0.479308067841920, 0.260345966079040, 0.}, 0.087807628716604},
  {{0.065130102902216, 0.065130102902216, 0.}, 0.026673617804419},
  {{0.065130102902216, 0.869739794195568, 0.}, 0.026673617804419},
  {{0.869739794195568, 0.065130102902216, 0.}, 0.026673617804419},
  {{0.312865496004874, 0.638444188569810, 0.}, 0.038556880445128},
  {{0.638444188569810, 0.048690315425316, 0.}, 0.038556880445128},
  {{0.048690315425316, 0.312865496004874, 0.}, 0.038556880445128},
  {{0.312865496004874, 0.048690315425316, 0.}, 0.038556880445128},
  {{0.638444188569810, 0.312865496004874, 0.}, 0.038556880445128},
  {{0.048690315425316, 0.638444188569810, 0.}, 0.038556880445128}};
// 1 negative weight, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 8 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP8Solin[16] = {
  {{0.333333333333333, 0.333333333333333, 0.}, 0.072157803838894},
  {{0.459292588292723, 0.459292588292723, 0.}, 0.047545817133643},
  {{0.459292588292723, 0.081414823414554, 0.}, 0.047545817133643},
  {{0.081414823414554, 0.459292588292723, 0.}, 0.047545817133643},
  {{0.170569307751760, 0.170569307751760, 0.}, 0.051608685267359},
  {{0.170569307751760, 0.658861384496480, 0.}, 0.051608685267359},
  {{0.658861384496480, 0.170569307751760, 0.}, 0.051608685267359},
  {{0.050547228317031, 0.050547228317031, 0.}, 0.016229248811599},
  {{0.050547228317031, 0.898905543365938, 0.}, 0.016229248811599},
  {{0.898905543365938, 0.050547228317031, 0.}, 0.016229248811599},
  {{0.263112829634638, 0.728492392955404, 0.}, 0.013615157087217},
  {{0.728492392955404, 0.008394777409958, 0.}, 0.013615157087217},
  {{0.008394777409958, 0.263112829634638, 0.}, 0.013615157087217},
  {{0.263112829634638, 0.008394777409958, 0.}, 0.013615157087217},
  {{0.728492392955404, 0.263112829634638, 0.}, 0.013615157087217},
  {{0.008394777409958, 0.728492392955404, 0.}, 0.013615157087217}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 9 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP9Solin[19] = {
  {{0.3333333333330, 0.3333333333333, 0.}, 0.04856789814140},
  {{0.4896825191990, 0.4896825191990, 0.}, 0.01566735011355},
  {{0.4896825191990, 0.0206349616025, 0.}, 0.01566735011355},
  {{0.0206349616025, 0.4896825191990, 0.}, 0.01566735011355},
  {{0.4370895914930, 0.4370895914930, 0.}, 0.03891377050240},
  {{0.4370895914930, 0.1258208170140, 0.}, 0.03891377050240},
  {{0.1258208170140, 0.4370895914930, 0.}, 0.03891377050240},
  {{0.1882035356190, 0.1882035356190, 0.}, 0.03982386946360},
  {{0.1882035356190, 0.6235929287620, 0.}, 0.03982386946360},
  {{0.6235929287620, 0.1882035356190, 0.}, 0.03982386946360},
  {{0.0447295133945, 0.0447295133945, 0.}, 0.01278883782935},
  {{0.0447295133945, 0.9105409732110, 0.}, 0.01278883782935},
  {{0.9105409732110, 0.0447295133945, 0.}, 0.01278883782935},
  {{0.2219629891610, 0.7411985987840, 0.}, 0.02164176968865},
  {{0.7411985987840, 0.0368384120547, 0.}, 0.02164176968865},
  {{0.0368384120547, 0.2219629891610, 0.}, 0.02164176968865},
  {{0.2219629891610, 0.0368384120547, 0.}, 0.02164176968865},
  {{0.7411985987840, 0.2219629891610, 0.}, 0.02164176968865},
  {{0.0368384120547, 0.7411985987840, 0.}, 0.02164176968865}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 10 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP10Solin[25] = {
  {{0.3333333333330, 0.3333333333330, 0.}, 0.04540899519140},
  {{0.4855776333840, 0.4855776333840, 0.}, 0.01836297887825},
  {{0.4855776333840, 0.0288447332327, 0.}, 0.01836297887825},
  {{0.0288447332327, 0.4855776333840, 0.}, 0.01836297887825},
  {{0.1094815754850, 0.1094815754850, 0.}, 0.02266052971775},
  {{0.1094815754850, 0.7810368490300, 0.}, 0.02266052971775},
  {{0.7810368490300, 0.1094815754850, 0.}, 0.02266052971775},
  {{0.3079398387640, 0.5503529418210, 0.}, 0.03637895842270},
  {{0.5503529418210, 0.1417072194150, 0.}, 0.03637895842270},
  {{0.1417072194150, 0.3079398387640, 0.}, 0.03637895842270},
  {{0.3079398387640, 0.1417072194150, 0.}, 0.03637895842270},
  {{0.5503529418210, 0.3079398387640, 0.}, 0.03637895842270},
  {{0.1417072194150, 0.5503529418210, 0.}, 0.03637895842270},
  {{0.2466725606400, 0.7283239045970, 0.}, 0.01416362126555},
  {{0.7283239045970, 0.0250035347627, 0.}, 0.01416362126555},
  {{0.0250035347627, 0.2466725606400, 0.}, 0.01416362126555},
  {{0.2466725606400, 0.0250035347627, 0.}, 0.01416362126555},
  {{0.7283239045970, 0.2466725606400, 0.}, 0.01416362126555},
  {{0.0250035347627, 0.7283239045970, 0.}, 0.01416362126555},
  {{0.0668032510122, 0.9236559335870, 0.}, 0.00471083348185},
  {{0.9236559335870, 0.0095408154003, 0.}, 0.00471083348185},
  {{0.0095408154003, 0.0668032510122, 0.}, 0.00471083348185},
  {{0.0668032510122, 0.0095408154003, 0.}, 0.00471083348185},
  {{0.9236559335870, 0.0668032510122, 0.}, 0.00471083348185},
  {{0.0095408154003, 0.9236559335870, 0.}, 0.00471083348185}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 11 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP11Solin[27] = {
  {{+0.5346110482710, +0.5346110482710, 0.}, 0.00046350316448},
  {{-0.0692220965415, +0.5346110482710, 0.}, 0.00046350316448},
  {{+0.5346110482710, -0.0692220965415, 0.}, 0.00046350316448},
  {{+0.3989693029660, +0.3989693029660, 0.}, 0.03857476745740},
  {{+0.2020613940680, +0.3989693029660, 0.}, 0.03857476745740},
  {{+0.3989693029660, +0.2020613940680, 0.}, 0.03857476745740},
  {{+0.2033099004310, +0.2033099004310, 0.}, 0.02966148869040},
  {{+0.5933801991370, +0.2033099004310, 0.}, 0.02966148869040},
  {{+0.2033099004310, +0.5933801991370, 0.}, 0.02966148869040},
  {{+0.1193509122830, +0.1193509122830, 0.}, 0.01809227025170},
  {{+0.7612981754350, +0.1193509122830, 0.}, 0.01809227025170},
  {{+0.1193509122830, +0.7612981754350, 0.}, 0.01809227025170},
  {{+0.0323649481113, +0.0323649481113, 0.}, 0.00682986550135},
  {{+0.9352701037770, +0.0323649481113, 0.}, 0.00682986550135},
  {{+0.0323649481113, +0.9352701037770, 0.}, 0.00682986550135},
  {{+0.5932012134280, +0.3566206482610, 0.}, 0.02616855598110},
  {{+0.0501781383105, +0.5932012134280, 0.}, 0.02616855598110},
  {{+0.3566206482610, +0.0501781383105, 0.}, 0.02616855598110},
  {{+0.0501781383105, +0.3566206482610, 0.}, 0.02616855598110},
  {{+0.3566206482610, +0.5932012134280, 0.}, 0.02616855598110},
  {{+0.5932012134280, +0.0501781383105, 0.}, 0.02616855598110},
  {{+0.8074890031600, +0.1714889803040, 0.}, 0.01035382981955},
  {{+0.0210220165362, +0.8074890031600, 0.}, 0.01035382981955},
  {{+0.1714889803040, +0.0210220165362, 0.}, 0.01035382981955},
  {{+0.0210220165362, +0.1714889803040, 0.}, 0.01035382981955},
  {{+0.1714889803040, +0.8074890031600, 0.}, 0.01035382981955},
  {{+0.8074890031600, +0.0210220165362, 0.}, 0.01035382981955}};
// 0 negative weights, 3 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 12 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP12Solin[33] = {
  {{0.4882173897740, 0.4882173897740, 0.}, 0.01286553322025},
  {{0.4882173897740, 0.0235652204524, 0.}, 0.01286553322025},
  {{0.0235652204524, 0.4882173897740, 0.}, 0.01286553322025},
  {{0.4397243922940, 0.4397243922940, 0.}, 0.02184627226900},
  {{0.4397243922940, 0.1205512154110, 0.}, 0.02184627226900},
  {{0.1205512154110, 0.4397243922940, 0.}, 0.02184627226900},
  {{0.2712103850120, 0.2712103850120, 0.}, 0.03142911210895},
  {{0.2712103850120, 0.4575792299760, 0.}, 0.03142911210895},
  {{0.4575792299760, 0.2712103850120, 0.}, 0.03142911210895},
  {{0.1275761455420, 0.1275761455420, 0.}, 0.01739805646535},
  {{0.1275761455420, 0.7448477089170, 0.}, 0.01739805646535},
  {{0.7448477089170, 0.1275761455420, 0.}, 0.01739805646535},
  {{0.0213173504532, 0.0213173504532, 0.}, 0.00308313052578},
  {{0.0213173504532, 0.9573652990940, 0.}, 0.00308313052578},
  {{0.9573652990940, 0.0213173504532, 0.}, 0.00308313052578},
  {{0.2757132696860, 0.6089432357800, 0.}, 0.02018577888320},
  {{0.6089432357800, 0.1153434945350, 0.}, 0.02018577888320},
  {{0.1153434945350, 0.2757132696860, 0.}, 0.02018577888320},
  {{0.2757132696860, 0.1153434945350, 0.}, 0.02018577888320},
  {{0.6089432357800, 0.2757132696860, 0.}, 0.02018577888320},
  {{0.1153434945350, 0.6089432357800, 0.}, 0.02018577888320},
  {{0.2813255809900, 0.6958360867880, 0.}, 0.01117838660115},
  {{0.6958360867880, 0.0228383322223, 0.}, 0.01117838660115},
  {{0.0228383322223, 0.2813255809900, 0.}, 0.01117838660115},
  {{0.2813255809900, 0.0228383322223, 0.}, 0.01117838660115},
  {{0.6958360867880, 0.2813255809900, 0.}, 0.01117838660115},
  {{0.0228383322223, 0.6958360867880, 0.}, 0.01117838660115},
  {{0.1162519159080, 0.8580140335440, 0.}, 0.00865811555435},
  {{0.8580140335440, 0.0257340505483, 0.}, 0.00865811555435},
  {{0.0257340505483, 0.1162519159080, 0.}, 0.00865811555435},
  {{0.1162519159080, 0.0257340505483, 0.}, 0.00865811555435},
  {{0.8580140335440, 0.1162519159080, 0.}, 0.00865811555435},
  {{0.0257340505483, 0.8580140335440, 0.}, 0.00865811555435}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 13 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP13Solin[37] = {
  {{0.33333333333330, 0.33333333333330, 0.}, 0.026260461700400},
  {{0.49504818494000, 0.49504818494000, 0.}, 0.005640072604650},
  {{0.00990363012059, 0.49504818494000, 0.}, 0.005640072604650},
  {{0.49504818494000, 0.00990363012059, 0.}, 0.005640072604650},
  {{0.46871663511000, 0.46871663511000, 0.}, 0.015711759181250},
  {{0.06256672978090, 0.46871663511000, 0.}, 0.015711759181250},
  {{0.46871663511000, 0.06256672978090, 0.}, 0.015711759181250},
  {{0.41452133680100, 0.41452133680100, 0.}, 0.023536251252100},
  {{0.17095732639700, 0.41452133680100, 0.}, 0.023536251252100},
  {{0.41452133680100, 0.17095732639700, 0.}, 0.023536251252100},
  {{0.22939957204300, 0.22939957204300, 0.}, 0.023681793268200},
  {{0.54120085591400, 0.22939957204300, 0.}, 0.023681793268200},
  {{0.22939957204300, 0.54120085591400, 0.}, 0.023681793268200},
  {{0.11442449519600, 0.11442449519600, 0.}, 0.015583764522900},
  {{0.77115100960700, 0.11442449519600, 0.}, 0.015583764522900},
  {{0.11442449519600, 0.77115100960700, 0.}, 0.015583764522900},
  {{0.02481139136350, 0.02481139136350, 0.}, 0.003987885732535},
  {{0.95037721727300, 0.02481139136350, 0.}, 0.003987885732535},
  {{0.02481139136350, 0.95037721727300, 0.}, 0.003987885732535},
  {{0.63635117456200, 0.26879499705900, 0.}, 0.018424201364350},
  {{0.09485382837960, 0.63635117456200, 0.}, 0.018424201364350},
  {{0.26879499705900, 0.09485382837960, 0.}, 0.018424201364350},
  {{0.09485382837960, 0.26879499705900, 0.}, 0.018424201364350},
  {{0.26879499705900, 0.63635117456200, 0.}, 0.018424201364350},
  {{0.63635117456200, 0.09485382837960, 0.}, 0.018424201364350},
  {{0.69016915998700, 0.29173006673400, 0.}, 0.008700731651900},
  {{0.01810077327880, 0.69016915998700, 0.}, 0.008700731651900},
  {{0.29173006673400, 0.01810077327880, 0.}, 0.008700731651900},
  {{0.01810077327880, 0.29173006673400, 0.}, 0.008700731651900},
  {{0.29173006673400, 0.69016915998700, 0.}, 0.008700731651900},
  {{0.69016915998700, 0.01810077327880, 0.}, 0.008700731651900},
  {{0.85140953783400, 0.12635738549200, 0.}, 0.007760893419500},
  {{0.02223307667410, 0.85140953783400, 0.}, 0.007760893419500},
  {{0.12635738549200, 0.02223307667410, 0.}, 0.007760893419500},
  {{0.02223307667410, 0.12635738549200, 0.}, 0.007760893419500},
  {{0.12635738549200, 0.85140953783400, 0.}, 0.007760893419500},
  {{0.85140953783400, 0.02223307667410, 0.}, 0.007760893419500}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 14 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP14Solin[42] = {
  {{0.48896391036200, 0.48896391036200, 0.}, 0.01094179068470},
  {{0.02207217927560, 0.48896391036200, 0.}, 0.01094179068470},
  {{0.48896391036200, 0.02207217927560, 0.}, 0.01094179068470},
  {{0.41764471934000, 0.41764471934000, 0.}, 0.01639417677205},
  {{0.16471056131900, 0.41764471934000, 0.}, 0.01639417677205},
  {{0.41764471934000, 0.16471056131900, 0.}, 0.01639417677205},
  {{0.27347752830900, 0.27347752830900, 0.}, 0.02588705225365},
  {{0.45304494338200, 0.27347752830900, 0.}, 0.02588705225365},
  {{0.27347752830900, 0.45304494338200, 0.}, 0.02588705225365},
  {{0.17720553241300, 0.17720553241300, 0.}, 0.02108129436850},
  {{0.64558893517500, 0.17720553241300, 0.}, 0.02108129436850},
  {{0.17720553241300, 0.64558893517500, 0.}, 0.02108129436850},
  {{0.06179988309090, 0.06179988309090, 0.}, 0.00721684983490},
  {{0.87640023381800, 0.06179988309090, 0.}, 0.00721684983490},
  {{0.06179988309090, 0.87640023381800, 0.}, 0.00721684983490},
  {{0.01939096124870, 0.01939096124870, 0.}, 0.00246170180120},
  {{0.96121807750300, 0.01939096124870, 0.}, 0.00246170180120},
  {{0.01939096124870, 0.96121807750300, 0.}, 0.00246170180120},
  {{0.77060855477500, 0.17226668782100, 0.}, 0.01233287660630},
  {{0.05712475740360, 0.77060855477500, 0.}, 0.01233287660630},
  {{0.17226668782100, 0.05712475740360, 0.}, 0.01233287660630},
  {{0.05712475740360, 0.17226668782100, 0.}, 0.01233287660630},
  {{0.17226668782100, 0.77060855477500, 0.}, 0.01233287660630},
  {{0.77060855477500, 0.05712475740360, 0.}, 0.01233287660630},
  {{0.57022229084700, 0.33686145979600, 0.}, 0.01928575539355},
  {{0.09291624935700, 0.57022229084700, 0.}, 0.01928575539355},
  {{0.33686145979600, 0.09291624935700, 0.}, 0.01928575539355},
  {{0.09291624935700, 0.33686145979600, 0.}, 0.01928575539355},
  {{0.33686145979600, 0.57022229084700, 0.}, 0.01928575539355},
  {{0.57022229084700, 0.09291624935700, 0.}, 0.01928575539355},
  {{0.68698016780800, 0.29837288213600, 0.}, 0.00721815405675},
  {{0.01464695005570, 0.68698016780800, 0.}, 0.00721815405675},
  {{0.29837288213600, 0.01464695005570, 0.}, 0.00721815405675},
  {{0.01464695005570, 0.29837288213600, 0.}, 0.00721815405675},
  {{0.29837288213600, 0.68698016780800, 0.}, 0.00721815405675},
  {{0.68698016780800, 0.01464695005570, 0.}, 0.00721815405675},
  {{0.87975717137000, 0.11897449769700, 0.}, 0.00250511441925},
  {{0.00126833093287, 0.87975717137000, 0.}, 0.00250511441925},
  {{0.11897449769700, 0.00126833093287, 0.}, 0.00250511441925},
  {{0.00126833093287, 0.11897449769700, 0.}, 0.00250511441925},
  {{0.11897449769700, 0.87975717137000, 0.}, 0.00250511441925},
  {{0.87975717137000, 0.00126833093287, 0.}, 0.00250511441925}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 15 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP15Solin[48] = {
  {{+0.5069729168580, +0.5069729168580, 0.}, 0.000958437821425},
  {{-0.0139458337165, +0.5069729168580, 0.}, 0.000958437821425},
  {{+0.5069729168580, -0.0139458337165, 0.}, 0.000958437821425},
  {{+0.4314063542830, +0.4314063542830, 0.}, 0.022124513635550},
  {{+0.1371872914340, +0.4314063542830, 0.}, 0.022124513635550},
  {{+0.4314063542830, +0.1371872914340, 0.}, 0.022124513635550},
  {{+0.2776936448470, +0.2776936448470, 0.}, 0.025593274359450},
  {{+0.4446127103060, +0.2776936448470, 0.}, 0.025593274359450},
  {{+0.2776936448470, +0.4446127103060, 0.}, 0.025593274359450},
  {{+0.1264648910410, +0.1264648910410, 0.}, 0.011843867935350},
  {{+0.7470702179170, +0.1264648910410, 0.}, 0.011843867935350},
  {{+0.1264648910410, +0.7470702179170, 0.}, 0.011843867935350},
  {{+0.0708083859747, +0.0708083859747, 0.}, 0.006644887845000},
  {{+0.8583832280510, +0.0708083859747, 0.}, 0.006644887845000},
  {{+0.0708083859747, +0.8583832280510, 0.}, 0.006644887845000},
  {{+0.0189651702411, +0.0189651702411, 0.}, 0.002374458304095},
  {{+0.9620696595180, +0.0189651702411, 0.}, 0.002374458304095},
  {{+0.0189651702411, +0.9620696595180, 0.}, 0.002374458304095},
  {{+0.6049544668930, +0.2613113711400, 0.}, 0.019275036299800},
  {{+0.1337341619670, +0.6049544668930, 0.}, 0.019275036299800},
  {{+0.2613113711400, +0.1337341619670, 0.}, 0.019275036299800},
  {{+0.1337341619670, +0.2613113711400, 0.}, 0.019275036299800},
  {{+0.2613113711400, +0.6049544668930, 0.}, 0.019275036299800},
  {{+0.6049544668930, +0.1337341619670, 0.}, 0.019275036299800},
  {{+0.5755865555130, +0.3880467670900, 0.}, 0.013607907160300},
  {{+0.0363666773969, +0.5755865555130, 0.}, 0.013607907160300},
  {{+0.3880467670900, +0.0363666773969, 0.}, 0.013607907160300},
  {{+0.0363666773969, +0.3880467670900, 0.}, 0.013607907160300},
  {{+0.3880467670900, +0.5755865555130, 0.}, 0.013607907160300},
  {{+0.5755865555130, +0.0363666773969, 0.}, 0.013607907160300},
  {{+0.7244626630770, +0.2857122200500, 0.}, 0.001091038683400},
  {{-0.0101748831266, +0.7244626630770, 0.}, 0.001091038683400},
  {{+0.2857122200500, -0.0101748831266, 0.}, 0.001091038683400},
  {{-0.0101748831266, +0.2857122200500, 0.}, 0.001091038683400},
  {{+0.2857122200500, +0.7244626630770, 0.}, 0.001091038683400},
  {{+0.7244626630770, -0.0101748831266, 0.}, 0.001091038683400},
  {{+0.7475564660520, +0.2155996640720, 0.}, 0.010752659923850},
  {{+0.0368438698759, +0.7475564660520, 0.}, 0.010752659923850},
  {{+0.2155996640720, +0.0368438698759, 0.}, 0.010752659923850},
  {{+0.0368438698759, +0.2155996640720, 0.}, 0.010752659923850},
  {{+0.2155996640720, +0.7475564660520, 0.}, 0.010752659923850},
  {{+0.7475564660520, +0.0368438698759, 0.}, 0.010752659923850},
  {{+0.8839645740920, +0.1035756165760, 0.}, 0.003836971315525},
  {{+0.0124598093312, +0.8839645740920, 0.}, 0.003836971315525},
  {{+0.1035756165760, +0.0124598093312, 0.}, 0.003836971315525},
  {{+0.0124598093312, +0.1035756165760, 0.}, 0.003836971315525},
  {{+0.1035756165760, +0.8839645740920, 0.}, 0.003836971315525},
  {{+0.8839645740920, +0.0124598093312, 0.}, 0.003836971315525}};
// 0 negative weights, 9 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 16 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP16Solin[52] = {
  {{+0.33333333333330, +0.33333333333330, 0.}, 0.023437848713800},
  {{+0.49738054194800, +0.49738054194800, 0.}, 0.003202939289290},
  {{+0.00523891610312, +0.49738054194800, 0.}, 0.003202939289290},
  {{+0.49738054194800, +0.00523891610312, 0.}, 0.003202939289290},
  {{+0.41346943854900, +0.41346943854900, 0.}, 0.020855148369700},
  {{+0.17306112290100, +0.41346943854900, 0.}, 0.020855148369700},
  {{+0.41346943854900, +0.17306112290100, 0.}, 0.020855148369700},
  {{+0.47045859906700, +0.47045859906700, 0.}, 0.013445742125050},
  {{+0.05908280186600, +0.47045859906700, 0.}, 0.013445742125050},
  {{+0.47045859906700, +0.05908280186600, 0.}, 0.013445742125050},
  {{+0.24055374997000, +0.24055374997000, 0.}, 0.021066261380800},
  {{+0.51889250006100, +0.24055374997000, 0.}, 0.021066261380800},
  {{+0.24055374997000, +0.51889250006100, 0.}, 0.021066261380800},
  {{+0.14796579422300, +0.14796579422300, 0.}, 0.015000133421400},
  {{+0.70406841155500, +0.14796579422300, 0.}, 0.015000133421400},
  {{+0.14796579422300, +0.70406841155500, 0.}, 0.015000133421400},
  {{+0.07546518765750, +0.07546518765750, 0.}, 0.007100049462500},
  {{+0.84906962468500, +0.07546518765750, 0.}, 0.007100049462500},
  {{+0.07546518765750, +0.84906962468500, 0.}, 0.007100049462500},
  {{+0.01659640262300, +0.01659640262300, 0.}, 0.001791231175635},
  {{+0.96680719475400, +0.01659640262300, 0.}, 0.001791231175635},
  {{+0.01659640262300, +0.96680719475400, 0.}, 0.001791231175635},
  {{+0.59986871117500, +0.29655559658000, 0.}, 0.016386573730300},
  {{+0.10357569224500, +0.59986871117500, 0.}, 0.016386573730300},
  {{+0.29655559658000, +0.10357569224500, 0.}, 0.016386573730300},
  {{+0.10357569224500, +0.29655559658000, 0.}, 0.016386573730300},
  {{+0.29655559658000, +0.59986871117500, 0.}, 0.016386573730300},
  {{+0.59986871117500, +0.10357569224500, 0.}, 0.016386573730300},
  {{+0.64219352494200, +0.33772306340300, 0.}, 0.007649153124200},
  {{+0.02008341165540, +0.64219352494200, 0.}, 0.007649153124200},
  {{+0.33772306340300, +0.02008341165540, 0.}, 0.007649153124200},
  {{+0.02008341165540, +0.33772306340300, 0.}, 0.007649153124200},
  {{+0.33772306340300, +0.64219352494200, 0.}, 0.007649153124200},
  {{+0.64219352494200, +0.02008341165540, 0.}, 0.007649153124200},
  {{+0.79959272097100, +0.20474828164300, 0.}, 0.001193122096420},
  {{-0.00434100261414, +0.79959272097100, 0.}, 0.001193122096420},
  {{+0.20474828164300, -0.00434100261414, 0.}, 0.001193122096420},
  {{-0.00434100261414, +0.20474828164300, 0.}, 0.001193122096420},
  {{+0.20474828164300, +0.79959272097100, 0.}, 0.001193122096420},
  {{+0.79959272097100, -0.00434100261414, 0.}, 0.001193122096420},
  {{+0.76869972140100, +0.18935849213100, 0.}, 0.009542396377950},
  {{+0.04194178646800, +0.76869972140100, 0.}, 0.009542396377950},
  {{+0.18935849213100, +0.04194178646800, 0.}, 0.009542396377950},
  {{+0.04194178646800, +0.18935849213100, 0.}, 0.009542396377950},
  {{+0.18935849213100, +0.76869972140100, 0.}, 0.009542396377950},
  {{+0.76869972140100, +0.04194178646800, 0.}, 0.009542396377950},
  {{+0.90039906408700, +0.08528361568270, 0.}, 0.003425027273270},
  {{+0.01431732023070, +0.90039906408700, 0.}, 0.003425027273270},
  {{+0.08528361568270, +0.01431732023070, 0.}, 0.003425027273270},
  {{+0.01431732023070, +0.08528361568270, 0.}, 0.003425027273270},
  {{+0.08528361568270, +0.90039906408700, 0.}, 0.003425027273270},
  {{+0.90039906408700, +0.01431732023070, 0.}, 0.003425027273270}};
// 0 negative weights, 6 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 17 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP17Solin[61] = {
  {{0.33333333333330, 0.33333333333330, 0.}, 0.016718599645400},
  {{0.49717054055700, 0.49717054055700, 0.}, 0.002546707720255},
  {{0.00565891888645, 0.49717054055700, 0.}, 0.002546707720255},
  {{0.49717054055700, 0.00565891888645, 0.}, 0.002546707720255},
  {{0.48217632262500, 0.48217632262500, 0.}, 0.007335432263800},
  {{0.03564735475080, 0.48217632262500, 0.}, 0.007335432263800},
  {{0.48217632262500, 0.03564735475080, 0.}, 0.007335432263800},
  {{0.45023996902100, 0.45023996902100, 0.}, 0.012175439176850},
  {{0.09952006195840, 0.45023996902100, 0.}, 0.012175439176850},
  {{0.45023996902100, 0.09952006195840, 0.}, 0.012175439176850},
  {{0.40026623937700, 0.40026623937700, 0.}, 0.015553775434500},
  {{0.19946752124500, 0.40026623937700, 0.}, 0.015553775434500},
  {{0.40026623937700, 0.19946752124500, 0.}, 0.015553775434500},
  {{0.25214126797100, 0.25214126797100, 0.}, 0.015628555609300},
  {{0.49571746405800, 0.25214126797100, 0.}, 0.015628555609300},
  {{0.25214126797100, 0.49571746405800, 0.}, 0.015628555609300},
  {{0.16204700465800, 0.16204700465800, 0.}, 0.012407827169850},
  {{0.67590599068300, 0.16204700465800, 0.}, 0.012407827169850},
  {{0.16204700465800, 0.67590599068300, 0.}, 0.012407827169850},
  {{0.07587588226070, 0.07587588226070, 0.}, 0.007028036535300},
  {{0.84824823547900, 0.07587588226070, 0.}, 0.007028036535300},
  {{0.07587588226070, 0.84824823547900, 0.}, 0.007028036535300},
  {{0.01565472696780, 0.01565472696780, 0.}, 0.001597338086890},
  {{0.96869054606400, 0.01565472696780, 0.}, 0.001597338086890},
  {{0.01565472696780, 0.96869054606400, 0.}, 0.001597338086890},
  {{0.65549320380900, 0.33431986736400, 0.}, 0.004059827659495},
  {{0.01018692882690, 0.65549320380900, 0.}, 0.004059827659495},
  {{0.33431986736400, 0.01018692882690, 0.}, 0.004059827659495},
  {{0.01018692882690, 0.33431986736400, 0.}, 0.004059827659495},
  {{0.33431986736400, 0.65549320380900, 0.}, 0.004059827659495},
  {{0.65549320380900, 0.01018692882690, 0.}, 0.004059827659495},
  {{0.57233759053200, 0.29222153779700, 0.}, 0.013402871141600},
  {{0.13544087167100, 0.57233759053200, 0.}, 0.013402871141600},
  {{0.29222153779700, 0.13544087167100, 0.}, 0.013402871141600},
  {{0.13544087167100, 0.29222153779700, 0.}, 0.013402871141600},
  {{0.29222153779700, 0.57233759053200, 0.}, 0.013402871141600},
  {{0.57233759053200, 0.13544087167100, 0.}, 0.013402871141600},
  {{0.62600119028600, 0.31957488542300, 0.}, 0.009229996605400},
  {{0.05442392429060, 0.62600119028600, 0.}, 0.009229996605400},
  {{0.31957488542300, 0.05442392429060, 0.}, 0.009229996605400},
  {{0.05442392429060, 0.31957488542300, 0.}, 0.009229996605400},
  {{0.31957488542300, 0.62600119028600, 0.}, 0.009229996605400},
  {{0.62600119028600, 0.05442392429060, 0.}, 0.009229996605400},
  {{0.79642721497400, 0.19070422419200, 0.}, 0.004238434267165},
  {{0.01286856083360, 0.79642721497400, 0.}, 0.004238434267165},
  {{0.19070422419200, 0.01286856083360, 0.}, 0.004238434267165},
  {{0.01286856083360, 0.19070422419200, 0.}, 0.004238434267165},
  {{0.19070422419200, 0.79642721497400, 0.}, 0.004238434267165},
  {{0.79642721497400, 0.01286856083360, 0.}, 0.004238434267165},
  {{0.75235100593800, 0.18048321164900, 0.}, 0.009146398385000},
  {{0.06716578241350, 0.75235100593800, 0.}, 0.009146398385000},
  {{0.18048321164900, 0.06716578241350, 0.}, 0.009146398385000},
  {{0.06716578241350, 0.18048321164900, 0.}, 0.009146398385000},
  {{0.18048321164900, 0.75235100593800, 0.}, 0.009146398385000},
  {{0.75235100593800, 0.06716578241350, 0.}, 0.009146398385000},
  {{0.90462550409600, 0.08071131367960, 0.}, 0.003332816002085},
  {{0.01466318222480, 0.90462550409600, 0.}, 0.003332816002085},
  {{0.08071131367960, 0.01466318222480, 0.}, 0.003332816002085},
  {{0.01466318222480, 0.08071131367960, 0.}, 0.003332816002085},
  {{0.08071131367960, 0.90462550409600, 0.}, 0.003332816002085},
  {{0.90462550409600, 0.01466318222480, 0.}, 0.003332816002085}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 18 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP18Solin[70] = {
  {{+0.33333333333330, +0.33333333333330, 0.}, +0.015404969968800},
  {{+0.49334480863100, +0.49334480863100, 0.}, +0.004536218339700},
  {{+0.01331038273820, +0.49334480863100, 0.}, +0.004536218339700},
  {{+0.49334480863100, +0.01331038273820, 0.}, +0.004536218339700},
  {{+0.46921059424200, +0.46921059424200, 0.}, +0.009380658469800},
  {{+0.06157881151610, +0.46921059424200, 0.}, +0.009380658469800},
  {{+0.46921059424200, +0.06157881151610, 0.}, +0.009380658469800},
  {{+0.43628139588700, +0.43628139588700, 0.}, +0.009720548992750},
  {{+0.12743720822600, +0.43628139588700, 0.}, +0.009720548992750},
  {{+0.43628139588700, +0.12743720822600, 0.}, +0.009720548992750},
  {{+0.39484617067300, +0.39484617067300, 0.}, +0.013876974305400},
  {{+0.21030765865300, +0.39484617067300, 0.}, +0.013876974305400},
  {{+0.39484617067300, +0.21030765865300, 0.}, +0.013876974305400},
  {{+0.24979456880300, +0.24979456880300, 0.}, +0.016128112675750},
  {{+0.50041086239400, +0.24979456880300, 0.}, +0.016128112675750},
  {{+0.24979456880300, +0.50041086239400, 0.}, +0.016128112675750},
  {{+0.16143219374400, +0.16143219374400, 0.}, +0.012537016308450},
  {{+0.67713561251200, +0.16143219374400, 0.}, +0.012537016308450},
  {{+0.16143219374400, +0.67713561251200, 0.}, +0.012537016308450},
  {{+0.07659822748540, +0.07659822748540, 0.}, +0.007635963985900},
  {{+0.84680354502900, +0.07659822748540, 0.}, +0.007635963985900},
  {{+0.07659822748540, +0.84680354502900, 0.}, +0.007635963985900},
  {{+0.02425243935350, +0.02425243935350, 0.}, +0.003396961011480},
  {{+0.95149512129300, +0.02425243935350, 0.}, +0.003396961011480},
  {{+0.02425243935350, +0.95149512129300, 0.}, +0.003396961011480},
  {{+0.04314636721700, +0.04314636721700, 0.}, -0.001111549364960},
  {{+0.91370726556600, +0.04314636721700, 0.}, -0.001111549364960},
  {{+0.04314636721700, +0.91370726556600, 0.}, -0.001111549364960},
  {{+0.63265796885700, +0.35891149494100, 0.}, +0.003165957038205},
  {{+0.00843053620242, +0.63265796885700, 0.}, +0.003165957038205},
  {{+0.35891149494100, +0.00843053620242, 0.}, +0.003165957038205},
  {{+0.00843053620242, +0.35891149494100, 0.}, +0.003165957038205},
  {{+0.35891149494100, +0.63265796885700, 0.}, +0.003165957038205},
  {{+0.63265796885700, +0.00843053620242, 0.}, +0.003165957038205},
  {{+0.57441097151100, +0.29440247675200, 0.}, +0.013628769024550},
  {{+0.13118655173700, +0.57441097151100, 0.}, +0.013628769024550},
  {{+0.29440247675200, +0.13118655173700, 0.}, +0.013628769024550},
  {{+0.13118655173700, +0.29440247675200, 0.}, +0.013628769024550},
  {{+0.29440247675200, +0.57441097151100, 0.}, +0.013628769024550},
  {{+0.57441097151100, +0.13118655173700, 0.}, +0.013628769024550},
  {{+0.62477904679300, +0.32501780164200, 0.}, +0.008838392824750},
  {{+0.05020315156570, +0.62477904679300, 0.}, +0.008838392824750},
  {{+0.32501780164200, +0.05020315156570, 0.}, +0.008838392824750},
  {{+0.05020315156570, +0.32501780164200, 0.}, +0.008838392824750},
  {{+0.32501780164200, +0.62477904679300, 0.}, +0.008838392824750},
  {{+0.62477904679300, +0.05020315156570, 0.}, +0.008838392824750},
  {{+0.74893317652300, +0.18473755966600, 0.}, +0.009189742319050},
  {{+0.06632926381090, +0.74893317652300, 0.}, +0.009189742319050},
  {{+0.18473755966600, +0.06632926381090, 0.}, +0.009189742319050},
  {{+0.06632926381090, +0.18473755966600, 0.}, +0.009189742319050},
  {{+0.18473755966600, +0.74893317652300, 0.}, +0.009189742319050},
  {{+0.74893317652300, +0.06632926381090, 0.}, +0.009189742319050},
  {{+0.76920700542000, +0.21879680001300, 0.}, +0.004052366404095},
  {{+0.01199619456620, +0.76920700542000, 0.}, +0.004052366404095},
  {{+0.21879680001300, +0.01199619456620, 0.}, +0.004052366404095},
  {{+0.01199619456620, +0.21879680001300, 0.}, +0.004052366404095},
  {{+0.21879680001300, +0.76920700542000, 0.}, +0.004052366404095},
  {{+0.76920700542000, +0.01199619456620, 0.}, +0.004052366404095},
  {{+0.88396230227300, +0.10117959713600, 0.}, +0.003817064535365},
  {{+0.01485810059010, +0.88396230227300, 0.}, +0.003817064535365},
  {{+0.10117959713600, +0.01485810059010, 0.}, +0.003817064535365},
  {{+0.01485810059010, +0.10117959713600, 0.}, +0.003817064535365},
  {{+0.10117959713600, +0.88396230227300, 0.}, +0.003817064535365},
  {{+0.88396230227300, +0.01485810059010, 0.}, +0.003817064535365},
  {{+1.01434726001000, +0.02087475528260, 0.}, +2.3093830397e-05},
  {{-0.03522201528790, +1.01434726001000, 0.}, +2.3093830397e-05},
  {{+0.02087475528260, -0.03522201528790, 0.}, +2.3093830397e-05},
  {{-0.03522201528790, +0.02087475528260, 0.}, +2.3093830397e-05},
  {{+0.02087475528260, +1.01434726001000, 0.}, +2.3093830397e-05},
  {{+1.01434726001000, -0.03522201528790, 0.}, +2.3093830397e-05}};
// 3 negative weights, 6 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 19 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP19Solin[73] = {
  {{0.33333333333330, 0.33333333333330, 0.}, 0.016453165694450},
  {{0.48960998707300, 0.48960998707300, 0.}, 0.005165365945650},
  {{0.02078002585400, 0.48960998707300, 0.}, 0.005165365945650},
  {{0.48960998707300, 0.02078002585400, 0.}, 0.005165365945650},
  {{0.45453689269800, 0.45453689269800, 0.}, 0.011193623631500},
  {{0.09092621460420, 0.45453689269800, 0.}, 0.011193623631500},
  {{0.45453689269800, 0.09092621460420, 0.}, 0.011193623631500},
  {{0.40141668064900, 0.40141668064900, 0.}, 0.015133062934750},
  {{0.19716663870100, 0.40141668064900, 0.}, 0.015133062934750},
  {{0.40141668064900, 0.19716663870100, 0.}, 0.015133062934750},
  {{0.25555165440300, 0.25555165440300, 0.}, 0.015245483901100},
  {{0.48889669119400, 0.25555165440300, 0.}, 0.015245483901100},
  {{0.25555165440300, 0.48889669119400, 0.}, 0.015245483901100},
  {{0.17707794215200, 0.17707794215200, 0.}, 0.012079606370800},
  {{0.64584411569600, 0.17707794215200, 0.}, 0.012079606370800},
  {{0.17707794215200, 0.64584411569600, 0.}, 0.012079606370800},
  {{0.11006105322800, 0.11006105322800, 0.}, 0.008025401793400},
  {{0.77987789354400, 0.11006105322800, 0.}, 0.008025401793400},
  {{0.11006105322800, 0.77987789354400, 0.}, 0.008025401793400},
  {{0.05552862425180, 0.05552862425180, 0.}, 0.004042290130890},
  {{0.88894275149600, 0.05552862425180, 0.}, 0.004042290130890},
  {{0.05552862425180, 0.88894275149600, 0.}, 0.004042290130890},
  {{0.01262186377720, 0.01262186377720, 0.}, 0.001039681013740},
  {{0.97475627244600, 0.01262186377720, 0.}, 0.001039681013740},
  {{0.01262186377720, 0.97475627244600, 0.}, 0.001039681013740},
  {{0.60063379479500, 0.39575478735700, 0.}, 0.001942438452490},
  {{0.00361141784841, 0.60063379479500, 0.}, 0.001942438452490},
  {{0.39575478735700, 0.00361141784841, 0.}, 0.001942438452490},
  {{0.00361141784841, 0.39575478735700, 0.}, 0.001942438452490},
  {{0.39575478735700, 0.60063379479500, 0.}, 0.001942438452490},
  {{0.60063379479500, 0.00361141784841, 0.}, 0.001942438452490},
  {{0.55760326158900, 0.30792998388000, 0.}, 0.012787080306000},
  {{0.13446675453100, 0.55760326158900, 0.}, 0.012787080306000},
  {{0.30792998388000, 0.13446675453100, 0.}, 0.012787080306000},
  {{0.13446675453100, 0.30792998388000, 0.}, 0.012787080306000},
  {{0.30792998388000, 0.55760326158900, 0.}, 0.012787080306000},
  {{0.55760326158900, 0.13446675453100, 0.}, 0.012787080306000},
  {{0.72098702581700, 0.26456694840700, 0.}, 0.004440451786670},
  {{0.01444602577610, 0.72098702581700, 0.}, 0.004440451786670},
  {{0.26456694840700, 0.01444602577610, 0.}, 0.004440451786670},
  {{0.01444602577610, 0.26456694840700, 0.}, 0.004440451786670},
  {{0.26456694840700, 0.72098702581700, 0.}, 0.004440451786670},
  {{0.72098702581700, 0.01444602577610, 0.}, 0.004440451786670},
  {{0.59452706895600, 0.35853935220600, 0.}, 0.008062273380850},
  {{0.04693357883820, 0.59452706895600, 0.}, 0.008062273380850},
  {{0.35853935220600, 0.04693357883820, 0.}, 0.008062273380850},
  {{0.04693357883820, 0.35853935220600, 0.}, 0.008062273380850},
  {{0.35853935220600, 0.59452706895600, 0.}, 0.008062273380850},
  {{0.59452706895600, 0.04693357883820, 0.}, 0.008062273380850},
  {{0.83933147368100, 0.15780740596900, 0.}, 0.001245970908745},
  {{0.00286112035057, 0.83933147368100, 0.}, 0.001245970908745},
  {{0.15780740596900, 0.00286112035057, 0.}, 0.001245970908745},
  {{0.00286112035057, 0.15780740596900, 0.}, 0.001245970908745},
  {{0.15780740596900, 0.83933147368100, 0.}, 0.001245970908745},
  {{0.83933147368100, 0.00286112035057, 0.}, 0.001245970908745},
  {{0.70108797892600, 0.07505059697590, 0.}, 0.009121420059500},
  {{0.22386142409800, 0.70108797892600, 0.}, 0.009121420059500},
  {{0.07505059697590, 0.22386142409800, 0.}, 0.009121420059500},
  {{0.22386142409800, 0.07505059697590, 0.}, 0.009121420059500},
  {{0.07505059697590, 0.70108797892600, 0.}, 0.009121420059500},
  {{0.70108797892600, 0.22386142409800, 0.}, 0.009121420059500},
  {{0.82293132407000, 0.14242160111300, 0.}, 0.005129281868100},
  {{0.03464707481680, 0.82293132407000, 0.}, 0.005129281868100},
  {{0.14242160111300, 0.03464707481680, 0.}, 0.005129281868100},
  {{0.03464707481680, 0.14242160111300, 0.}, 0.005129281868100},
  {{0.14242160111300, 0.82293132407000, 0.}, 0.005129281868100},
  {{0.82293132407000, 0.03464707481680, 0.}, 0.005129281868100},
  {{0.92434425262100, 0.06549462808290, 0.}, 0.001899964427650},
  {{0.01016111929630, 0.92434425262100, 0.}, 0.001899964427650},
  {{0.06549462808290, 0.01016111929630, 0.}, 0.001899964427650},
  {{0.01016111929630, 0.06549462808290, 0.}, 0.001899964427650},
  {{0.06549462808290, 0.92434425262100, 0.}, 0.001899964427650},
  {{0.92434425262100, 0.01016111929630, 0.}, 0.001899964427650}};
// 0 negative weights, 0 points outside of the triangle,  total sum of the
// weights is 0.5

// -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 20 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */

IntPt triP20Solin[79] = {
  {{+0.3333333333333, +0.3333333333333, 0.}, +0.016528527770800},
  {{+0.5009504643520, +0.5009504643520, 0.}, +0.000433509592831},
  {{-0.0019009287044, +0.5009504643520, 0.}, +0.000433509592831},
  {{+0.5009504643520, -0.0019009287044, 0.}, +0.000433509592831},
  {{+0.4882129579350, +0.4882129579350, 0.}, +0.005830026358200},
  {{+0.0235740841305, +0.4882129579350, 0.}, +0.005830026358200},
  {{+0.4882129579350, +0.0235740841305, 0.}, +0.005830026358200},
  {{+0.4551366869500, +0.4551366869500, 0.}, +0.011438468178200},
  {{+0.0897266360994, +0.4551366869500, 0.}, +0.011438468178200},
  {{+0.4551366869500, +0.0897266360994, 0.}, +0.011438468178200},
  {{+0.4019962593180, +0.4019962593180, 0.}, +0.015224491336950},
  {{+0.1960074813630, +0.4019962593180, 0.}, +0.015224491336950},
  {{+0.4019962593180, +0.1960074813630, 0.}, +0.015224491336950},
  {{+0.2558929097590, +0.2558929097590, 0.}, +0.015312445862700},
  {{+0.4882141804810, +0.2558929097590, 0.}, +0.015312445862700},
  {{+0.2558929097590, +0.4882141804810, 0.}, +0.015312445862700},
  {{+0.1764882559950, +0.1764882559950, 0.}, +0.012184028838400},
  {{+0.6470234880100, +0.1764882559950, 0.}, +0.012184028838400},
  {{+0.1764882559950, +0.6470234880100, 0.}, +0.012184028838400},
  {{+0.1041708553370, +0.1041708553370, 0.}, +0.007998716016000},
  {{+0.7916582893260, +0.1041708553370, 0.}, +0.007998716016000},
  {{+0.1041708553370, +0.7916582893260, 0.}, +0.007998716016000},
  {{+0.0530689638409, +0.0530689638409, 0.}, +0.003849150907800},
  {{+0.8938620723180, +0.0530689638409, 0.}, +0.003849150907800},
  {{+0.0530689638409, +0.8938620723180, 0.}, +0.003849150907800},
  {{+0.0416187151960, +0.0416187151960, 0.}, -0.000316030248744},
  {{+0.9167625696080, +0.0416187151960, 0.}, -0.000316030248744},
  {{+0.0416187151960, +0.9167625696080, 0.}, -0.000316030248744},
  {{+0.0115819214068, +0.0115819214068, 0.}, +0.000875567150595},
  {{+0.9768361571860, +0.0115819214068, 0.}, +0.000875567150595},
  {{+0.0115819214068, +0.9768361571860, 0.}, +0.000875567150595},
  {{+0.6064026461060, +0.3448557702290, 0.}, +0.008232919594800},
  {{+0.0487415836648, +0.6064026461060, 0.}, +0.008232919594800},
  {{+0.3448557702290, +0.0487415836648, 0.}, +0.008232919594800},
  {{+0.0487415836648, +0.3448557702290, 0.}, +0.008232919594800},
  {{+0.3448557702290, +0.6064026461060, 0.}, +0.008232919594800},
  {{+0.6064026461060, +0.0487415836648, 0.}, +0.008232919594800},
  {{+0.6158426144570, +0.3778432695950, 0.}, +0.002419516770245},
  {{+0.0063141159486, +0.6158426144570, 0.}, +0.002419516770245},
  {{+0.3778432695950, +0.0063141159486, 0.}, +0.002419516770245},
  {{+0.0063141159486, +0.3778432695950, 0.}, +0.002419516770245},
  {{+0.3778432695950, +0.6158426144570, 0.}, +0.002419516770245},
  {{+0.6158426144570, +0.0063141159486, 0.}, +0.002419516770245},
  {{+0.5590480003900, +0.3066354790620, 0.}, +0.012902453267350},
  {{+0.1343165205470, +0.5590480003900, 0.}, +0.012902453267350},
  {{+0.3066354790620, +0.1343165205470, 0.}, +0.012902453267350},
  {{+0.1343165205470, +0.3066354790620, 0.}, +0.012902453267350},
  {{+0.3066354790620, +0.5590480003900, 0.}, +0.012902453267350},
  {{+0.5590480003900, +0.1343165205470, 0.}, +0.012902453267350},
  {{+0.7366067432630, +0.2494193627750, 0.}, +0.004235545527220},
  {{+0.0139738939624, +0.7366067432630, 0.}, +0.004235545527220},
  {{+0.2494193627750, +0.0139738939624, 0.}, +0.004235545527220},
  {{+0.0139738939624, +0.2494193627750, 0.}, +0.004235545527220},
  {{+0.2494193627750, +0.7366067432630, 0.}, +0.004235545527220},
  {{+0.7366067432630, +0.0139738939624, 0.}, +0.004235545527220},
  {{+0.7116751422870, +0.2127757248030, 0.}, +0.009177457053150},
  {{+0.0755491329098, +0.7116751422870, 0.}, +0.009177457053150},
  {{+0.2127757248030, +0.0755491329098, 0.}, +0.009177457053150},
  {{+0.0755491329098, +0.2127757248030, 0.}, +0.009177457053150},
  {{+0.2127757248030, +0.7116751422870, 0.}, +0.009177457053150},
  {{+0.7116751422870, +0.0755491329098, 0.}, +0.009177457053150},
  {{+0.8614027171550, +0.1469654360530, 0.}, +0.000352202338954},
  {{-0.0083681532082, +0.8614027171550, 0.}, +0.000352202338954},
  {{+0.1469654360530, -0.0083681532082, 0.}, +0.000352202338954},
  {{-0.0083681532082, +0.1469654360530, 0.}, +0.000352202338954},
  {{+0.1469654360530, +0.8614027171550, 0.}, +0.000352202338954},
  {{+0.8614027171550, -0.0083681532082, 0.}, +0.000352202338954},
  {{+0.8355869579120, +0.1377269788290, 0.}, +0.005056342463750},
  {{+0.0266860632587, +0.8355869579120, 0.}, +0.005056342463750},
  {{+0.1377269788290, +0.0266860632587, 0.}, +0.005056342463750},
  {{+0.0266860632587, +0.1377269788290, 0.}, +0.005056342463750},
  {{+0.1377269788290, +0.8355869579120, 0.}, +0.005056342463750},
  {{+0.8355869579120, +0.0266860632587, 0.}, +0.005056342463750},
  {{+0.9297561715570, +0.0596961091490, 0.}, +0.001786954692975},
  {{+0.0105477192941, +0.9297561715570, 0.}, +0.001786954692975},
  {{+0.0596961091490, +0.0105477192941, 0.}, +0.001786954692975},
  {{+0.0105477192941, +0.0596961091490, 0.}, +0.001786954692975},
  {{+0.0596961091490, +0.9297561715570, 0.}, +0.001786954692975},
  {{+0.9297561715570, +0.0105477192941, 0.}, +0.001786954692975}};
// 3 negative weights, 9 points outside of the triangle,  total sum of the
// weights is 0.5

static IntPt *GQTSolin[21] = {
  triP1Solin,  triP1Solin,  triP2Solin,  triP3Solin,  triP4Solin,  triP5Solin,
  triP6Solin,  triP7Solin,  triP8Solin,  triP9Solin,  triP10Solin, triP11Solin,
  triP12Solin, triP13Solin, triP14Solin, triP15Solin, triP16Solin, triP17Solin,
  triP18Solin, triP19Solin, triP20Solin};
static int GQTnPtSolin[21] = {1,  1,  3,  4,  6,  7,  12, 13, 16, 19, 25,
                              27, 33, 37, 42, 48, 52, 61, 70, 73, 79};
static std::vector<IntPt *> GQTGL(40, nullptr);

IntPt *getGQTPts(int order)
{
  if(order < 21) return GQTSolin[order];
  int n = (order + 3) / 2;
  if(static_cast<int>(GQTGL.size()) < order + 1)
    GQTGL.resize(order + 1, nullptr);
  if(!GQTGL[order]) {
    int npts = n * n;
    IntPt *intpt = new IntPt[npts];
    GaussLegendreTri(n, n, intpt);
    GQTGL[order] = intpt;
  }
  return GQTGL[order];
}

int getNGQTPts(int order)
{
  if(order < 21) return GQTnPtSolin[order];
  return ((order + 3) / 2) * ((order + 3) / 2);
}
